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Hydrodynamic interactions in a suspension of spherical particles confined between 
two parallel planar walls are studied under creeping-flow conditions. The many-particle 
friction matrix in this system is evaluated using our novel numerical algorithm based 
on transformations between Cartesian and spherical representations of Stokes flow. The 
Cartesian representation is used to describe the interaction of the fluid with the walls 
and the spherical representation is used to describe the interaction with the particles. 
The transformations between these two representations are given in a closed form, which 
allows us to evaluate the coefficients in linear equations for the induced-force multipoles 
on particle surfaces. The friction matrix is obtained from these equations, supplemented 
with the superposition lubrication corrections. We have used our algorithm to evalu- 
ate the friction matrix for a single sphere, a pair of spheres, and for linear chains of 
spheres. The friction matrix exhibits a crossover from a quasi-two-dimensional behavior 
(for systems with small wall separation H) to the three-dimensional behavior (when the 
distance H is much larger than the interparticle distance L). The crossover is especially 
pronounced for a long chain moving in the direction normal to its orientation and parallel 
to the walls. In this configuration, a large pressure buildup occurs in front of the chain 
for small values of the gapwidth _ff , which results in a large hydrodynamic friction force. 
A standard wall superposition approximation does not capture this behavior. 



1. Introduction 

Numerous recent papers reflect a growing interest in the static and dynamic proper- 
ties of suspensions in confined geometries. There are investigations o f the formation of 
colloidal crystals on pat terned and planar surfaces fLin et aL"20nnt ISeehg et a/.l 120021: 
ISubram anian e t al\\l99^ . studies of single- file diffusion of Brownian particles in a chan- 
nel jWei et aLll200Cl|) . an d experiments on quasi- two-dimensional suspensions confined 
between two planar walls JCarb aial-Tin qco et a/.lll997tlLancon et ai]l200ll:lMarcus et all 
0^£; Santana-Solano fc Arauz-Lara 200ll). Quasi- two-dimensional suspensions of parti- 
cles adsorbed at a fluid interface llZahn et oill997tlR,inn et QiJll999tlCichocki et oli2004) 
or confined in a thin liquid film ijSethumadhavan. Nikolov fc Wasanll200lj) have also been 
examined. 

Experiments on quasi-two-dimensional systems revealed many striking phenomena 

like, for i nstance, the first-order transitions between fluid, hexatic, and solid phases 

I Marcus fc Ric6il997.1 . string-hke cooperative motion of suspension particles ijMarcus. Schofield fc Ricel 

t On leave from IPPT Warsaw, Poland 
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1^22 )j Eind oscillatory melting of a crystalline phase in shear flow (fStancik fc Hawkins 

200.1^ . Other interesting examples include a logarithmic behavior of the mean-s quare dis- 

)lacen ient of Brownian pa rticles in quasi - two-di mensional systems, predicted bv lCichocki fc Felderhoj 
1994|) and observed by iMarcus et al\ (Il999tl: a hy drodynamic enhancement of self- 



migration of partic les in Poiseuille flow towards the channel center ijNott fc Brad" 



Q) , and 

2uM 



diffusion for strongly-charged particles ijzahn et a/.l 1l997: .Pesc hc Nagclc 200Q,), and 
migration of partic l 
iLvon fc Leailll99j^ . 

The particle-wall and interparticle interaction potentials fully determine the equilib- 
rium structure of confined colloidal suspensions. The dynamics of such systems, how- 
ever, is significantly affected by the many-body hydrodynamic forces. For spherical par- 



tides in an unbounded space, efficient algorithms for evaluat 
mobility matrices have been developed llDurlofskv. Brady 


ing many-body frictio 
& Bossi.cJ 119871 iT.addl 


n and 

mi 


ICichocki. Felderhof. Hinsen. Wainrvb & BlawzdziewicJl994 


:ISieron & Rradv'2nni 


).Us- 



ing the image representation technique, such algorithms have been generalized for parti- 
cles adsorbed on a planar fl uid-air interface llCichocki et al. 12004) and for particles con- 



fined in a thin liquid film (jBlawzdziewicz fcWainrvb[]200^ ). The image-representation 



method has also been proposed for a suspen sion bounded by a single rigid planar wall 
" Cichocki fc Joneslll998HCichocki et a^Jl2000|) . A ;wo-wall generalization of this method 



Bhattacharya fc Blawzdzi ewicd 12002 oj) and several other techniques were used to de- 



scribe motion of an individual p article between two planar walls ( Ganatos et al 19806IqI : 

IStaben et Qi1l2003H.Tonesll2004l) . 

The hydrodynamics of many particles in the two- wall geometry is much more complex, 
and available results are limited. Durlofsk v fc Bradv. (,198 9) have developed a method 
that combines boundary-integral and Stokesian-dynamics elements. In their approach, 
the walls are discretized, and the particles are represented using low-order force mul- 
tipoles and lubrication contributions, as in the standard Stokesian-dynamics algorithm 
I Durlofskv et oi.l^l987^ It seems that this m etho d has not been further pursued. In an 



alternative approach, Nott fc Bradvl l|l994l) and iMorris fc Bradvl l)l998^ studied flows 



in wall-bounded suspensions by modeling the walls as static, closely packed arrays of 
spheres, and using the standard Stokesian-dynamics algorithm for an unbounded system 
to evaluate the motion of the suspended particles. The results obtained in this way are 
only qualitative, especially for small wall separations, because the wal ls are porous and 
rough . Recently, a two-wall superposition approximation was used by IPesche fc Nagelel 
l|2000l) and several other groups, but the validity range of this approximation cannot be 
determined without comparison with more accurate results. 

In our paper we present a novel, highly accurate algorithm to evaluate the many-body 
hydrodynamic interactions in a suspension of spherical particles confined between two 
planar walls. In our approach, the flow fleld in the system is expanded using two basis 
sets of solutions of Stokes equations — the spherical and Cartesian bases. The spherical 
basis is applied to describe the flow field scattered from the particles, and the Cartesian 
basis is used in the analysis of the interaction of the fiow with the walls. The key result 
of our study is a set of transformation formulas for conversion between the spherical 
and Cartesian representations. In our algorithm the expansion of the fiow field into the 
basis fields is combined with the two-particle superposition approximation for the friction 
matrix in order to incorporate slowly convergent lubrication corrections. Since the force 
multipoles induced on particle surfaces are included to arbitrary order, results of arbitrary 
accuracy are obtained. 

Our paper is organized as follows. The induced-force formulation of the problem is de- 
scribed in a-nd the multipolar representation of the flow in terms of force multipoles 
induced on the particles is recalled in ^ Our main theoretical results are outlined in ^ 
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and ^ The Cartesian basis set of Stokes flows is defined in 21 along with the transfor- 
mation relations between the Cartesian and spherical bases, displacement theorems for 
the Cartesian basis fields and expressions for the wall-reflection matrix. These essential 
elements are combined in ^to evaluate the wall contribution to the Green's matrix. The 
numerical implementation of our method is outlined in ^ Examples of numerical results 
(for a single particle, two particles and many-particle systems) are provided in SjTl The 
multiparticle results have been selected to illustrate the crossover between the quasi-two- 
dimensional and three-dimensional behavior of the friction matrix as a function of the 
interparticle distance. 

Since the full description of the theory underlying our algorithm requires more space, 
this paper outlines only the most important elements and lists the crucial results indis- 
pensable for the numerical implementation. The details of our theoretical analysis and 
a more complete description of the algorithm are presented in a separate publication 
l)Bhattacharva et a/...2005.'l , hereafter referred to as Ref. I. 



2. Multiparticle hydro dynamic interactions 

2.1. Hydrodynamic resistance 

We consider a suspension of N spherical particles of the radius a in a creeping flow 
between two parallel planar walls. The no-slip boundary conditions are assumed on the 
particles and on the walls. The walls are at the positions z = and z — H, where H is the 
separation between walls, and r = (x, y, z) are the Cartesian coordinates. The position of 
the center of particle i is denoted by — {Xi, Yi, Zi), the translational and rotational 
particle velocities are denoted by Ui and O^, and the external forces and torques acting 
on the particle are denoted by and Ti. 

We focus on a system of spheres undergoing translational and rotational rigid-body 
motion with no external flow. As in an unbounded space, the particle dynamics in the 
system is characterized by the resistance matrix 



1, 



defined by the linear relation 



N 

E 



(2.1) 



(2.2) 



between the translational and rotational particle velocities and the forces and torques. 
The dot in equation (|2.2(l denotes the matrix multiplication and contraction of the Carte- 
sian tensorial components of the resistance matrix. 



2.2. Induced-force formulation 

The effect of the suspended particles on the surrounding fluid can be described in terms 
of the induced force distributions on the particle surfaces 



where 



F,(r) =a-2<5(r,-a)f,(r). 



Hi 



(2.3) 
(2.4) 
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and ri = |ri|. By definition of the induced force, the flow field 

N 

v(r)=5] /T(r,r')-F,(r')dr' (2.5) 

2—1 

is identical to the velocity field in the presenc e of the moving particles ijCox fc BrenneJ 
llflfiTt iMaznr Redeanixlll 974 lFelderhollll973^ . Here 

T(r,r')=To(r-r') + T'(r,r') (2.6) 

is the Green's function for the Stokes flow in the presence of the boundaries; the Green's 
function T(r,r') is decomposed into the Oseen tensor To(r — r') and the part T'(r,r') 
that describes the flow reflected from the walls. In equation (|2.5|) it is assumed that the 
particles move with given velocities, but no external flow is imposed. 

The resistance relation (|2.2|l is linked to the induced-force distributions 1)2. 3|l through 
the expressions 

J^, = yF,(r)dr, T, = yr,:XF,(r)dr (2.7) 

for the total force and torque, respectively. To determine the resistance matrix 1)2. l|l 
we thus need to evaluate the induced forces (|2.3|) for given translational and angular 
velocities of the particles. 

2.3. Boundary-integral equations for the induced forces 

For a system of particles moving with the translational and angular velocities Ui and fl^, 
the induced-force distribution 1)2. 3|l can be obtained from the boundary-integral equation 
of the form 

N 

[Z-iF,](r)+^ /[(l-%)To(r-r') + T'(r,r')]-F,(r')dr'=vf(r), r e 5„ (2.8) 
where 

vf(r) X r, (2.9) 

is the rigid-body velocity fleld associated with the particle motion, and Si is the surface of 
particle i. In the boundary-integral equation H2.8|l . denotes the one-particle scattering 
operator that describes the response of an individual particle to an external flow in an 
unbounded space. This operator is defined by the linear relation 

F,--Z,(vr-vf), (2.10) 

where vj" is the velocity incident to part icle i. For specific particl e models, explicit 
expressions for the operat or Z^ are known ijjones fc Schmitz,.198& .Cichocki et aL.1988t 
iBlawzdziewicz et al\\l99^ . 



3. Force-multipole expansion 

3.1. Spherical basis fields 

As in a standard force-multipole approach ijCichocki et a/.lll994i l2000l) the boundary- 
integral equation ()2.8(l is transformed into a linear matrix equation by projecting it onto 
a spherical basis of S tokes flow. To this end we use the reciprocal basis sets defined by 
iGichocki et (Hlll988^ : we introduce, however, a slightly different normalization to exploit 
the full symmetry of the problem. 
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The singular and nonsingular spherical basis solutions of Stokes equations Vj~ (r) and 
'^tmai^) (with I — 1,2,...; m — —I, . . . , Z; and tr = 0, 1, 2) have the following separable 
form in the spherical coordinates r = (r, 0, 0): 

vr™.w = vr_(0,0K('+-), (3.1a) 



v,|„.(r)=V+„,(0,0)r'+-\ (3.1&) 

where the coefficients Vj~j^^(6', (j)) and Vj^^(6', (f>) are combinations of vector spherical har- 
monics with angular order I and azimuthal order m. This property and the r-dependence 
in equations H3.1|l define the Stokes-flow fields Vj^j^(r) up to a normalization constant. 
Explicit expressions for the functions Vj^^^^ in our present normalization are given in 
Appendix The iustificat i on for this choice of the normalization is discussed in Ref. L 
Following ICichocki et all ()l988(l we also introduce the reciprocal basis fields W;'^^(r), 
defined here by the orthogonality relations of the form 

^ Inia I ^ V ra' a' 



{€^La I ^Um'a') = <5.., , (3.2) 



where 

Sl{r)^a-'6{r-a), (3.3) 

and 

(A I B) = J A*(r) • B(r) dr. (3.4) 

The asterisk in the above relation denotes the complex conjugate. We note that due to 
the proper choice of defini ng properties of the spherical basis sets, the basis fields v^^^ 
and satisfy relation (jCichocki et a/.lll983) 

v,™.(r) = V J To(r - r')<5,S(r')w+,.(r') dr', (3.5) 

where r] is the viscosity of the fluid. Relation 1)3. 5|) assures that the Lorentz reciprocal 
symmetry of Stokes flow is reflected in th e symmetry of the resulting matrix representa- 
tion of the problem ijCichocki et aZ.ll200Cl|) . 

3.2. Matrix representation 

The matrix representation of the boundary- integral equation (|2.8II is obtained by expand- 
ing the force distributions induced on each particle into the induced- force multipoles. The 
force multipolar moments of the force distribution 1)2. 3|l are defined by the relation 

F,(r) = m^^^^'^Kn - a)w+_(r.), (3-6) 

Iraa 

where is the relative position 1)2.4)1 with respect to the particle center. According to 
equations ()3.5)) and 1)3. 6)1 . the multipolar moments fi{lma) are identical (apart from the 
trivial factor rj) to the expansion coefficients of the flow field scattered by an isolated 
particle in unbounded space into the singular basis fields v^^^. 

To obtain a linear matrix equation for the set of force multipolar moments fi{lma), 
the multipolar representation of the induced force ()3.6)l is inserted into the boundary- 
integral equation 1)2.8)1 . and the resulting expression is expanded into the nonsingular 
basis solutions ()3.16)l . In particular, for the rigid-body velocity field we have 

vf(r)=^c.(/™a)v+_(r,), (3.7) 
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where the expansion coefScients Ci{lma) are nonzero only for I — 1 and cr = 0, 1. 

As the resuh of this procedure we get the hnear force-multipole equation, which can 
be written in the form 

N 

Mijilm I l'm')-fj{l'm') = Ci{lm). (3.8) 

j — l I'm' 

We use here a matrix notation in the three-dimensional linear space with the compo- 
nents corresponding to the indices cr = 0, 1, 2 that identify the tensorial character of the 
basis flow fields (|3.1|) . Accordingly, the arrays fj{l'm') and Ci{lm) have the components 
fj{l'm'a') and Ci{lma) and the matrix My (?m | I'm') has the elements Mij{lma \ I'm! a'), 
where cr, cr' = 0, 1, 2. The many-particle resistance matrix (|2.1(l can be obtained by solv- 
ing equation 1)3. 8|l and projecting the induced force multipoles onto the total force and 
torque 1)2. 7|l . Explicit expressions for the resistance matrix in terms of the generalized 
friction matrix F = M^^ are given in Appendix IbI 

For a wall bounded system the matrix M can be decomposed into three contributions 

M„(Zm I I'm') = %<5,P^wZ-\0 + Q%{lm \ I'm') + Q',,{lm \ I'm'). (3.9) 

The first term Z~"^(?) corresponds to the one particle operator in equation (|2.8() . 
Accordingly, the matrix Zi(/) relates the force multipoles fi{l'm') induced on particle i 
to the coefficients in the expansion of the flow field incoming to this particle into the 
nonsingular spherical basis fields l|3.16|) . By spherical symmetry, this ter m is diagonal in 
the m ultipolar orders I and to, and for rigid spheres it is explicitly known ijCichocki et all 
119881^ . 

The Green matrices G^j(Zto | I'm') and G'ij{lm \ I'm') correspond to the integral 
operators with the k ernels T n(r — r') and T'(r, r ') in equation (|2.8(l . As discussed by 
ICichocki et al and bv iBhattacharva et all l|200!4l the matrix G^j{lm \ I'm') coin- 

cides (apart from the normalization f actors) with the displace ment matrix for spherical 
basis fields, which is explicitly known llFelderhof fc ,Toneslll989t) . The only unknown com- 
ponent in expression H3.9|l is thus the wall contribution 

G',^{lma I I'm'a') = (rOw+„^(r,) | V,„,,,,{v,)), (3.10) 

where 

v^,„,,,(r) = J T'(r,r')5f (r')w+„,,,(r')dr'. (3.11) 

In the following sections we express this contribution in terms of the Cartesian basis set 
of Stokes flows. 



4. Cartesian representation 

The difficulties associated with the evaluation of the matrix G'^ in systems with two 



planar walls stem from the incompatibility of the spherical basis fields with the wall 



geometry. In particular, the image representation of a force multipole I Cichocki fc Jone^ 
[l99 8; Bhattacharya & Blawzdziewicz 2002_a) is insufficient for a two-wall system, due 
to the slow convergence of the multiple-image series. We propose here an alternative 
technique, which relies on a transformation between the spherical basis fields (|3.1|l and 
a Cartesian basis set of Stokes flows. In the Cartesian representation the flow reflected 
from the walls can be obtained in a closed form; thus the difficulties associated with the 
multiple-image series are avoided. 

According to our Cartesian representation method, the wall contribution H3.10|l to 
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the matrix M is evaluated by (j) expanding the spherical basis flow field V;7„j,g./ (rj) 
produced by a force multipole at the position rj into the Cartesian basis; (ii) solving the 
two-wall problem in the Cartesian representation; and (Hi) transforming the resulting 
reflected flow back to the spherical basis (j3.1&|l centered at the position r^. As a result 
of this procedure, the matrix elements (|3.10() are expressed in terms of two-dimensional 
Fourier integrals with respect to the lateral coordinates x,y. Our method is outlined in 
the following sections. 

4.1. Cartesian basis 

To describe the flow field between two walls parallel to the x-y plane, it is convenient to 
use a basis set of Stokes flows of the separable form 

vJ.(r)=V±,(z)e"'-''±'=^ (4.1) 

that is consistent with the wall geometry. Here 

p = xe^ + ysy (4.2) 

is the projection of the vector r onto the x-y plane, 

k — ' ^) 

is the corresponding two-dimensional wave vector, and k = |k|. By analogy to the spher- 
ical basis (|3.1|l . there exist three types of solutions cr = 0, 1, 2 for each k. These solutions 
involve a potential flow, a flow with nonzero vorticity, and a pressure-driven flow. Explicit 
expressions for the coefficients V^^(z) are given in Appendix ICl To achieve a substantial 
simpliflcation of our final results, the relative amplitudes of the three components in the 
basis fields (|4.1|) have been carefully chosen, as discussed in Ref. I. 

4.2. Transformation relations 

The transformation relations between the spherical and Cartesian sets of solutions of 
Stokes equations can be expressed by the formulas 

vr™.W= / dk'^v±^,(r)T±s-(k',/m;a'k), 

±z < 0, (4.4) 

v^.(r) - E v,+™'.'Wrs+c^(^'™',k;a' I <y). (4.5) 

I'm' a' 

As demonstrated in Ref. I, the transformation matrices TQg"(k, Im) and Tg(^(/m, k) have 
the factorized form 

Ttiilm, k) = (-i)"(27rfc)-i/2g-im^K(fc, /) • t+^f (/m), (4.6a) 

T±<7(k, Im) = i'"(27rfc)-i/2gi™vt±-(;^) . K(fc, I), {4.6b) 

where tp is the polar angle in the Fourier space, 

K{k,l;a\a')^S^„'k'+''~\ (4.7) 

and the matrices Tg(^(/r7i) and 'f^^{lm) are independent of the wave vector k. Due 
to the proper choice of the spherical and Cartesian fields, the transformation matrices 
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Tg(^(Zm) and T^g {Im) have a simple symmetric form 



sc 



a b c 
2a 2b 
4a 



(4.8a) 



' sc 



= fe]^ = (-1) 



l-\-m 



(4.86) 



b a 
-2a 


where the dagger denotes the Hermitian conjugate. The three independent scalar coefR- 
cients in equations H4.8|) are 



c 
-2b 
4a 



a = - m)\{l + m)\{2l + l)]-l/^ 



(4.9a) 



b = 2am/l, 



(4.9&) 



1{2P ~2l-l)-2m'^{l-2) 
' l{2l - 1) ' 



(4.9c) 



4.3. Cartesian displacement formulas 

In an analysis of the flow between the walls it is convenient to use the Cartesian basis 
fields centered at different positions (e.g., the particle or wall position). As shown in Ref. 
I, the fields H4.I|I centered at different points Ri and R2 are related by the displacement 
formula 

v^.(r2) =5]v±,,(ri)5±±(Ri2,k;a' | a), (4.10) 

a' 

where ri = r — Ri, Y2 = r — R2, and R12 = Ri — R2. Since the shift of the origin 
of the coordinate system preserves the behavior of the flow fields H4.II) at infinity, the 
superscripts in equation (|4.I0() are either all positive or all negative. The displacement 
matrices S^^(Ri2,k) can be factorized as follows. 



Sg±(Ri2,k) 



S±±(fcZl2)e'''•Pl^ 



where 



1 
1 

-2kZ I 



-kZ 



S++{kZ) = 



1 2kZ 
I 
I 



,fcZ 



and 



R-12 — Pi2 + ^1262. 



(4.11) 



(4.12) 



(4.13) 



4.4. Single-wall reflection matrix 

The Cartesian basis fields (|4.I|I are well suited for a description of the interaction of 
the flow with planar walls because, due to the translational invariance of the problem, 
the lateral Fourier modes with different wave vectors k do not couple. The effect of a 
single wall on the flow field in the system can be characterized in terms of the one-wall 
reflection matrix Z^. To define this quantity we consider Stokes flow in a fluid bounded 
by a single wall in the plane 

z = Zw. (4.14) 
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The fluid occupies either the halfspace z > (denoted by $1+) or z < Zw (denoted by 

rj-). 

The velocity field in the halfspace fi^ can be uniquely decomposed into the incoming 
and reflected flows 

v(r)=v^^(r)+vr(r). (4.15) 
The flow v™(r) is nonsingular in the halfspace il^, and the flow v™*(r) is nonsingular in 
the halfspace fi^. Thus these flows have the following expansions in the Cartesian basis, 

v;"(r)= / dk^c|,-(ka)v±^(rw), (4.16a) 

vr(r)= / dk5]cr(kcT)vJ,(rw). (4.166) 

Here 

rw = r-Rw (4.17) 

denotes the position of the point r relative to the wall, where ~ (Xw,l'w,^w) has 
arbitrary lateral coordinates and Y^. 

The single- wall scattering matrix relates the expansion coefficients of the incoming 
and refiected flows: 

cr(k) = -Z„.C(k), (4.18) 
where c°"*(k) and cJJ,^(k) denote the arrays of expansion coefficients in equations (|4.16|) . 
For a rigid wall with no-slip boundary conditions we have 



1 
1 
1 



(4.19) 



as shown in Ref. I. For planar interfaces with other boundary conditions (e.g., a surfact ant- 
covered fluid-fluid interface discussed bv lBlawzdziewicz. Cristini fc Loewenberell999|) the 
scattering matrix can also be obtained. 



5. Evaluation of the wall contribution to Green's matrix 

5.1. Single-wall system 

The transformation, displacement, and reflection matrices described in 21 can be used 
to construct the matrix G'ij for a suspension bounded by a single planar wall or by two 
planar walls. For a single wall, the matrix H3.10|l can be expressed by the two-dimensional 
Fourier integral 

G\,ilm\l'm')^ y dk*s(k;Z,w,Zw,)e"'-''- (5.1) 

with the integrand of the form 

§s(k; Z™, Zw,) = -rj^'T+^ilm, k) • S^^{kZ,^) • • S^^{kZ^,) ■ T±g-(k, I'm'), (5.2) 

where Zi-^ = Zi — Z^ and Z^j = Z^ — Zj are the vertical offsets between the points i 
and j and the wall. 

The physical interpretation of relation H5.2|l is straightforward. First, the spherical 
components of the flow produced by a multipolar force distribution at point j are trans- 
formed by the matrix Tq^ into the Cartesian basis. The Cartesian components of the 
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velocity field are propagated by the matrix Sq^(Rwj") to the wall, where they are scat- 
tered, as represented by the matrix Zw The reflected field is propagated by the matrix 
SQ^(Riw) to the point i, and, finally, the flow is transformed by the matrix Tg^f back 
into the spherical basis. 

Due to symmetry properties of the component matrices (cf., relations (|4.6() . (|4.8() . 
H4.12|l . and (|4.19|) ) the wall contribution to the Green's matrix 1)5. l|l satisfies the Lorentz 
symmetry 

G'ijilm I I'm') = G'^lil'm' \ Im). (5.3) 

We note that for the single-wall problem the Fourier integr al (15.111 can be exp l icitly per- 
formed, which yields the image-singularity result derived bv lCichocki fc .Toned lll998^. As 
discussed in ^0 both the Fourier representation 1)5. l|l and the result of ICichocki fc JonesI 
lfl99R^ are used in our algorithm to accelerate the convergence of the two-wall integrals 
by a subtraction of the single-wall contributions. 

5.2. Two-wall system 

The single- wall result presented above can be generalized to the flow between two parallel 
walls. We assume that the walls are in the planes 

z = Zl, z = Zu, (5.4) 

where 

< Z^. (5.5) 

The two-wall Green's matrix H3.10|l can be expressed in the form analogous to equations 
l|K?T|l and i-e., 



G',j(lm\l'm')^ j dk*(k;Z,L,Z,L,^Lu)e'''''''^ 



(5.6) 



^'(k;Z,L,ZjL,^Lu) =-ry"'"rsc(/TO,k).S,w(k)-ZTw(k)-Swj(k)- Tcs(k,Z'm'), (5.7) 
where the component matrices are given by 

T+<7(k, Im) 



7cs(k, Zm) 

Tc^(k, Im) 
7sc(?m,k) = [ T+c (?m,k) T++(/m,k) ] , 

" S++(ZL,k) 

Swj(k) 

Sc-(^u,k) 



(5.8a) 
(5.8&) 
(5.9a) 



and 



Sjw(k) 



ZTw(k) 







S++(Z,uk) 

Z-i S++(ZLuk) 
Sc-(ZuLk) Z^i 



(5.9&) 



(5.10) 
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The physical interpretation of equation H5.7|l is similar to the interpretation of relation 
l|5.2(l . except that two separate sets of expansion coefficients are now used for the flow 
field incoming to the lower and upper walls, which is reflected in the corresponding 
block structure of matrices (|5.8|) - (|5.1U|) . The matrix Tcs(k, Z'm') transforms the field 
produced by a force multipole at the position Rj into the Cartesian basis; the basis 
fields are used in the region Z < Zj and the basis fields Vj^^ in the region Z > Zj, 
consistent with relation (|4.4|) . The Cartesian fields are then translated to the positions 
of the respective walls by the matrix Swj(k). The matrix ZTw(k), defined by equation 
H5.1U|I . describes the interaction of the flow field with both walls. We note that this matrix 
involves the displacements matrices Sq~ (Zui^li) and Sj'''(ZLuk), which correspond to 
the propagation of the flow field between the walls in the multiple-reflection process. 
After reflection from the walls is completed, the matrix Siw(k) propagates the flow field 
to the target point i, and the matrix Tsc('™,k) transforms it back into the spherical 
representation. 

Due to the symmetries of the 3x3 transformation and displacement matrices, the 
corresponding symmetry relations 

Tcs(k,M = [7"sca™,k)]t, (5.11a) 



Sw^(k) = [S,w(k)]^ (5.116) 

ZTw(k) = [ZTw(k)]^ (5.11c) 

are satisfied by the matrices (|5.8|I - H5.10|I . Equation H5.7f) thus implies that the Green 
matrix H5.6|l satisfies the Lorentz symmetry (15.3(1 . 



6. Numerical implementation 

The evaluation of the resistance matrix ^.^-^ from relations given in Appendix ^ re- 
quires solving the linear algebraic equation H3.8|l for the array of induced-force multipo- 
lar moments in order to obtain the generalized friction coefficients Fij{lma \ I'm' a'). 
In expression 1(3. 9|l for the matrix M^- the single particle scattering matrix and 
the unbounde d-spac e Green's matrix are known explicitly ( Felderh of fc Jo nes 198!^ 
Icichocki et a/.lll988|) . The remaining term — the two- wall contribution G'ij — is evaluated 
numerically, using relations H5.6|l - H5.10|l along with our expressions for the Cartesian dis- 
placement matrices (|4.12|l . the transformation matrices (I4.8|l . and the single- wall scat- 
tering matrix H4.19|l . 

Taking into account the structure H4.6|l of the transformation matrices Tg^T and , 
the angular integral in equation (|5.1|) can be performed analytically. The integration 
yields the result in the form of a Hankel transform of the order m' — m. Accordingly, 
only a one-dimensional integral in equation (|5.6|l has to be performed numerically. The 
numerical integration is straightforward when the lateral separation between particles i 
and j is small compared to the wall separation. For large interparticle separations pij, 
however, the integration is more difficult due to the oscillatory behavior of the integrand. 

To avoid numerical integration of a highly oscillatory function, the Fourier amplitude 
in (|5.()|) is decomposed 

vE'(k) = ^'L(k) -H*u(k) -l-(5*(k) (6.1) 

into the superposition of the single- wall contributions ^'l and and the remaining 
part representing hydrodynamic interactions between the walls. From an analysis of 
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expression ()5.2|l we find that the magnitude of the single- waU Fourier amplitudes 5'L(k) 
and vE'u(k) for large k is 

^'^(k) - e-'^'^^', a = L,U, (6.2) 

where A^"^ is the vertical offset between the point i and the reflection of point j in the 
wall a. In contrast, the large-fc behavior of the wall-interaction part of Fourier amplitude 
(jnH) is 

<5^(k) - e~'='^'^ (6.3) 

where 

Xj=2H~\Z^j\> H. (6.4) 
The lengthscale A^- equals the vertical offset \Zi — Z"| between the target point i and 
the closer of the two second-order images of the source point j. Since (S5'(k) decays on 
the wave- vector scale set by the distance between the walls H > min(A|^'' , A^^-*), a 
smaller number of oscillations of the Fourier modes contributes to the integral after the 
single-wall terms have been subtracted. 

In our algorithm, the short-range function (5^'(k) is integrated numerically. The one- 
particle contributions \&L(k) and 5'u(k) ar e evaluated analyt i callv, using the exphcit 
image- representation expressions derived in ICichocki fc Jones! ((l998|) . In this way we 
avoid integrating a highly oscillatory function when the particles are close to a wall. 
The procedure can be further improved, either by subt racting several terms associated 
with h igher-order wall reflections of the source multipole jBhattacharva fc Blawzdziewicd 
or by using asymptotic formulas for the integrals H5.6|l . We have recently derived 
a complete set of such expressions, which will be presented in a separate publication. 

In order to improve convergence with the order /max of the multipoles inclu ded in the 
calcul ation we employ a standard technique, originally introduced by Durlofskv et al\ 
l)l987j) . Accordingly, the lubrication forces that cause a slow convergence of the results 
with Imax are included in the friction matrix using a superposition approximation. Both, 
the interparticle and particle- wall lubrica tion corrections are inc luded in this way. Follow- 
ing the implementation of this method bv lCichocki et oil ( 200n|) for a single wall problem, 
we represent the elements of resistance matrix in the form 

C,-C-r' + Cr'" + ^Cr (6-5) 

Here CiJ'''^ denotes the superposition of two-particle resistance matrices evaluated for 
isolated particle pairs in the unbounded space, and 

C-r"=^.. E CUi)- (6.6) 

Q=L,U 

is the superposition of one-particle contributions in the presence of individual walls. 
The one-particle contributions can be evaluated using a series expansion of r esistance 
coefhcients in inverse powers of particle- wall separation ("Cichocki fc ,Io nealll998l) m com- 
bination with the appropriate lubrication results (Kim fc Karrila 1991). The two-particle 
superposition contributions C^j are evaluated in a similar way. The convergence with 
the multipolar truncation order /max for the quantity A^^ is fast; some convergence tests 
are presented in Ref. I. 

In the present implementation of our method, the numerical cost scales as 0{N^) with 
the number of particles N, because the linear equation 1)3. 9|l is solved by inversion of 
the matrix M. However, the numerical efficiency of our algorithm can be substantially 
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Figure 1. Lateral and vertical components of the translational friction matrix H7.8|l for a 
single sphere between two parallel walls, versus dimensionless gap 17.411 . Center particle position 
h = (solid line); off-center position h — (dashed line). Dotted lines represent the 
lubrication results 17.511 . 

improved by applying fast-multipole or PPPM acceleration methods in combination with 
asymptotic expressions for the elements of the matrix M. 

7. Results 

In this section we present a set of numerical results for hydrodynamic interactions in 
systems of spherical particles confined between two parallel planar walls. Our goal is 
both, to illustrate typical behavior of the hydrodynamic friction matrix for particles in 
the confined region, and to demonstrate the capabilities of our numerical algorithm. The 
results for a single particle and for pairs of particles, shown in figuresmSl were obtained 
using the multipolar approximation with the truncation at the order l^^x = 12. This 
truncation is sufficient to obtain results with the accuracy better than the resolution of 
the plots, even for the smallest wall separation H considered. The multi-particle results 
in figures IMTTl were obtained using l^^x — 8. 

7.1. Single particle 
In figure n the lateral and vertical friction coefhcients 

C||=c^r = c^?^ a-cir- (7.1) 

are shown for a single particle at the center and off-center positions 

h^\H, (7.2a, h) 

where h is the distance of the particle from the lower wall. The results are normalized 
by the Stokes friction coefficient Co — Gnrja, 

qi = Cii/Co, a = a/Co, (7.3) 

and are plotted versus the normalized gap 

e^h/a^l (7.4) 

between the particle and the lower wall. As expected, for small values of the gap, the 
lateral and vertical resistance coefficients approach the asymptotic lubrication behavior 
(in figure n indicated by dotted lines). For h = the lubrication behavior is 

Cli =-iloge-t-C(i), d = 2e-i, (7.5a) 
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where the singular terms correspond to the superpos ition of tw o particle-wall lubrication 
regions (cf. lubrication expressions given bv .Cichocki fc Jones, IPOS') . For the off-center 
position h =^ there is only one lubrication region, thus 

Cii =-^loge + C(i), a^(-\ {7.5b) 

A comparison of the numerical results shown in figure ^ with the asymptotic behavior 
CTSl yields C(i) = 1.45 and C(i) = 1.49. 

We n ote that our one-particle results agree with the numerical calculation by Ganatos et all 
l|l980fcl fah and with our earlier results |B hattacharva fc Blawzdziewicjl2 002fc ) obtained 



by an image- representation method I Bhattacharva fc Blawzdziewicjl2002tJ) . 



7.2. Two particles 

Sample results for the translational components of the two-particle resistance matrix 

C^^ = CT^/Co, *,J = 1,2, (7.6) 

(where a,/3 = x,y,z) are presented in figures [21|H1 The results in figures ^re shown 
for horizontal particle configurations hi — h2 — h, where hi is the distance of particle i 
from the lower wall. As for a single sphere, we consider the center and off-center positions 
H7.2|) . The relative horizontal displacement of the particles is 

Pl2 ^ PuGx. (7.7) 

To emphasize the crossover between the three-dimensional behavior for pi2 ^ H and 
a quasi-two-dimensional behavior for pi2 3> i?, we discuss our results in terms of the 
dimensionless variables scaled by the distance between the walls H, 

P = Pi2/H, d = a/H, L = p — 2a. {7.8a, b,c) 

The resistance coefficients in figures |2JE1 are plotted versus the dimensionless separation 
between the particle surfaces L. 

Self -resistance coefficients 

Figures 121 and O illustrate the behavior of the diagonal components of the translational 
self- and mutual resistance matrices and Ci'^, respectively, and figure^lshows the off- 
diagonal elements ^ff and . The remaining coefficients of the two-particle translational 
resistance matrix either vanish or can be related to the above coefficients by symmetry. 

The results for the resistance coefficients ("f presented in figure |21 are scaled by the 
corresponding single-particle friction coefficients (|7.3|l . For small distances between the 
particle surfaces L <^ a the longitudinal resistance coefficient Cif is dominated by the 
0{L^^) interparticle lubrication friction; the lubrication behavior of the components 
transverse to the direction of the line connecting the particle centers is C,\i, Cii ~ logi. 

In the intermediate regime p « 1 the two-particle friction matrix undergoes a crossover 
to a quasi-two-dimensional far-field asymptotic behavior at large interparticle distances. 
A signature of the crossover is the kink seen in the plot of C,\\ for the particles at the 
center position h = ^H. For large interparticle separations p ^ 1, the lateral components 
of the self-friction matrix approach the one-particle asymptotic value as 

Cn « CTf = Cll + O(r'), P » 1- (7.9) 

This result stems from the far-field behavior of the disturbance velocity produced by the 
particles. For the lateral motion the far-field disturbance decays as 0{p^^), as shown in 
Ref. I. Since the contribution of the second particle to the self-components of the friction 
matrix involves the field scattered back to the first particle, the asymptotic behavior 




Figure 2. Diagonal components of translational self-resistance matrix 17.61 for pair particles 
between two planar walls, scaled by corresponding one-particle values, versus dimensionless 
distance (17. Sb ) between particle surfaces. Walls are in planes z = 0,H, and particles are on 
axis X. The left and right panels correspond to center and off-center particle configurations (as 
indicated). Dimensionless gap between the particles and the closer wall e = 0.02 (solid line); 
e — 0.1 (long-dashed); e — 1.0 (dash-dotted); e = 4 (short-dashed); e = 16 (dotted). 



1)7.91) is obtained. In contrast, the disturbance field corresponding to the vertical motion 
decays exponentially, which yields an exponential approach of the vertical component of 
the friction matrix (ff to the one-particle value (± . 

Mutual resistance coefficients 

An analogous reasoning applied to the mutual components of the friction matrix yields 
the asymptotic behavior 

Cfl « -Ci2^ = Ap-' + O(r'), P » 1, (7.10) 
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Figure 3. Same as figure |51 except that for mutual components of the resistance matrix. The 
lateral components i^ff and ^12 a-re scaled by the asymptotic result = Ap~'^ corresponding 
to relation 17.101 . and the vertical component is CJl^ unsealed. 



where the amplitude A > depends on the size of the particles and on their vertical 
positions in the gap. Note that the sign of the transverse resistance coefficient ^12 ^-t 
large interparticle distances is opposite to the sign of the corresponding coefficient in the 
unbounded space. 

The results for C,^^ and Cfl shown in figure 13 are scaled using expression (|7.1Q(I , with 
the amplitude A plotted in figure |Sl(discussed below). Since decays exponentially for 
large p, the results for this component are presented unsealed. Similar to the results in 
figure |5] for the self-resistance matrix, the diagonal components of the mutual resistance 
matrix have a lubrication singularity for particles in contact, and for p = 0(1) they 
exhibit a crossover to the asymptotic 0{p~^) far-field behavior in the regime p ^ 1. The 
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Figure 4. Cross components of self- and mutual resistance matrix I7.6|l for the off-center 
particle configuration h = \H. Dimensionless particle-wall gaps e corresponding to different 
lines are the same as in figure |21 
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Figure 5. Amplitude A of the 0{p~^) far field asymptotic behavior H7.1()^ of the mutual 
resistance matrix, versus particle size normalized by distance between walls a = aj H . Center 
particle configuration /ii = /12 = \H (solid line); off center configuration hi — h2 = -^H (dashed 
line). Circles represent the corresponding asymptotic values IT.llt . 
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p/Aa (,11 i,^^ (,11 C,i2 <,12 '■12 

1.01 15.08 2.31 3.26 -13.75 -0.51 -0.24 

1.1 3.35 1.98 2.95 -1.99 -0.117 0.132 

2.0 1.93 1.87 2.80 -0.368 0.200 0.117 

5.0 1.85 1.85 2.79 -0.064 0.063 0.00015 

Table 1. Normalization factors 17.131 for the configurations represented in figures |S| and |7| 



near-field and far-field region are most pronounced for large values of the particle-wall 
gap e (i.e., for a ^ 1) because of the lengthscale separation. 

Cross-terms 

The cross-elements of the self- and mutual friction matrix and (^f| are shown 
(unsealed) in figure 01 Since for the center particle position 17.2b ) these components 
vanish by symmetry, the results are presented only for the off-center configuration (|7.2b '). 
The nonzero values of the cross-resistance coefficients C,fi and C,fi arise indirectly, due to 
the asymmetry of the flow field scattered from the walls. Therefore, for L — Q there is no 
lubrication singularity. The cross-resistance coefficients involve vertical particle motion; 
thus, for large interparticle separations and decay exponentially. 

Amplitude of the Jar-field asymptotic behavior 

The behavior (|7.9|l and (|7.1U|) of the two-particle re sistance coefficients for p ^ 1 is 
consistent with the asymptotic expressions derived by iLiron fc MochonI l[l976j) for the 
far-field flow produced by Stokeslets oriented in the direction parallel and normal to the 
walls. Using the Liron-Mochon expression and applying the Stokes resistance formula 
to evaluate forces acting on small particles in the space between the walls yields the 
asymptotic behavior (|7.1U|I . with the amplitude given by 

4 = 9/ii(l-/ii)/i2(l -^2) + 0(a), (7.11) 
a 

where hi = hi/H. Figure [S] shows the dependence of the far-field amplitude A on the 
dimensionless particle size a; the limiting result (|7.11(l is indicated by circles. 

We emphasize that the far-field form of the disturbance flow produced by particles in 
a domain bounded by parallel walls, and the corresponding properties of the resistance 
matrix are important for understanding the macroscopic dynamics of suspensions in 
slit-pore geometries. A more detailed analysis of this problem will be given elsewhere. 

Skew configurations 

So far we have focused on horizontal particle configurations with both particles at the 
same distance from the walls. In figures EHHl '^6 consider skew configurations with the 
vertical positions 

hi^h, h2=H-h. (7.12) 

Figures El and [7| show the diagonal components Cf" and (,^2 of the self- and mutual 
resistance matrices, and figure |S1 presents the off-diagonal components and C,fi. The 
results are plotted versus the normalized particle- wall gap \7A\ . 

The diagonal resistance coefRcients in figures and 13 are scaled by the value 
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Figure 6. Diagonal components of translational self-resistance matrix for skew configurations 
l|7.12|l of a particle pair between two walls separated by distance H/2a = 2, versus normalized 
particle-wall gap 17.41 . The results are scaled by the value 17.1311 for the center position, which 
corresponds to e = 1. Lateral particle separation p/2a — 1.01 (solid lines); 1.1 (dashed); 2 
(dash-dotted); 5 (dotted). 

corresponding to the center configuration of the particle pair at a given lateral separa- 
tion p and wall-to-wall distance H. The resistance coefficients for the center configuration 
H7.13|l have been discussed above; the values of the normalization factors for the param- 
eter values represented in figures and [T] are listed in tabled 

For large lateral interparticle distances, the self-resistance coefficients, shown in figure 
El approach the corresponding one-particle results. We note that due to the fast asymp- 
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Figure 7. Same as figure |S1 except that for mutual friction coefficients. 

totic approach H7.9|l . the resuhs for p/2a = 5 essentially coincide with the one-particle 
values. For small particle-wall gaps, the lateral coefficients ^if and (^ff exhibit the log- 
arithmic lubrication singularity, and the vertical component Cf f has the 1/e singularity. 
The rapid variation of the longitudinal coefficient Cif in the regime e ~ 1 (center particle 
positions) for p/2a = 1.01 results from the strong lubrication interaction between the 
particles. The same remark applies to the mutual longitudinal coefficient C12 shown in 
figure 

According to the results in table ^ the mutual resistance coefficients approach zero 
for large interparticle distances. The far-field behavior is consistent with the asymptotic 



Spherical particles between planar walls 



21 



0.0 

-0.2 

C~xz 
11 -0.4 

-0.6 

-0.8 

0.0 0.25 0.5 0.75 1.0 



1.0 
0.75 
Cl2 0.5 
0.25 



0.0 

0.0 0.25 0.5 0.75 1.0 

e 

Figure 8. Same as figure El except that for self- and mutual cross friction coefficients, and 

that the results are unsealed. 

expression (|7.1()() for the lateral components (^ff and and the asymptotic exponential 
decay for the vertical component Cff . We note that for small and moderate interparticle 
distances there is no simple relation between the components ^ff and d'^; however 
Ci2 ~ ~Ci2 for p 1, in agreement with equation (|7.10l) . 

The off-diagonal components Cff and (^f|, shown unsealed in figure|Hl are exponentially 
small for p ^ 1. For both particles at the center of the space between the walls these 
components vanish by symmetry. 

7.3. Multi-particle systems 

In figures l^llll we present some results for hydrodynamic resistance functions of rigid 
linear arrays of N touching spheres. The spheres are positioned on a line parallel to the 
axis X at the center of the space between the walls, i.e., 

h^^^H, i = l,...,N. (7.14) 

The diagonal components of the translational resistance matrix of the array treated as 
a single rigid body, evaluated per one sphere, 

N 

Cr = (^^Co)"' E C*r"' c. = x,y,z, (7.15) 
are plotted in figure |H| versus the number of spheres N in the chain. The results for the 
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Figure 9. Resistance coefficients per particle 17.1511 of rigid linear arrays of touching spheres 
on a line parallel to axis x at the center position H7.14|l . scaled by corresponding one-particle 
values 17.3II . versus number of spheres A'^. Dimensionless gap between the particles and walls 
e — 0.02 (solid circles); e = 0.1 (open circles); e = 1.0 (solid squares); e — 4 (open squares); 
e = 16 (solid triangles); e = oo (open triangles). 



longitudinal and transverse components and are shown normalized by the lateral 
one-particle resistance coefficient C||; the vertical component is normalized by (^±. 
The results indicate that for large separations between the walls, when compared to the 
chain length, all three components of the resistance matrix Cg" decrease monotonically 
with N, and behave as 1/ log for 1<C TV < H/2a. We also find that Cc" - Cc - "^Cc" 
in this regime (,Blawzdziewicz et oi...20n5.) . 
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Figure 10. Resistance coefRcients 17.16al l representing the total forces 17. 17b ) on individual 
spheres in the chain of length A'^ — 20, scaled by corresponding one-particle values l|7.3^ . versus 
particle number in the chain. The configuration of the chain and the dimensionless particle-wall 
gaps e are the same as in figure El 



For moderate and small values of the wall-to-wall distance H, however, the behavior 
of each component Cc" of the chain resistance matrix is qualitatively different. The longi- 
tudinal component decreases monotonically with iV, which is similar to the behavior 
in the unbounded space, but the variation is smaller. The vertical component ^^f initially 
increases with N, and then saturates at a constant value that depends on the wall sepa- 
ration H. In contrast, for small H, the transverse component increases linearly with 
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Figure 11. Same as figure IT?71 except that for resistance coeflicients Ij7.166^ representing the 
total torques 17.17b 'l on individual spheres, and that the results are plotted unsealed. Only half 
of the chain is shown; the resistance coefficients are antisymmetric with respect to the chain 
center. 

N in the range of chain lengths shown. Additional numerical simulations for chains with 
the length up to N=100 (not presented) indicate that the resistance coefficients even- 
tually saturate for large N. We note that the standard w all superposition approximation 
entirely misses this behavior ijBhattacharva et a/.ll2005|) . 

A better insight into the mechanisms underlying the above-illustrated qualitative fea- 
tures of the resistance matrix can be gained from the set of more detailed results for a 
chain of the length = 20 plotted in figures [TUI and ITTI In these figures, we show the 
resistance coefficients 

N 

Cr«=Co"'E^S""' c. = ^,y,^, (7.16a) 

and 

N 

C^"(z) = |(aCo)-^5^Ci;''", f3a = zy,yz, (7.166) 
i=i 

representing the normalized applied force and torque 

:F^ = e„Cr(*): = e^C^"(«) (7.17a, b) 

acting on particle i in a chain moving in the direction a with a unit velocity. By symmetry, 
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the forces act only in the direction of the chain motion, and the only nonzero torque 
coefficients are those listed in equation (|7.16&|l . 

According to the results shown in figure lTUI for the motion in the x direction, the forces 
acting on the first and the last particle in the chain are larger than the forces acting on 
the particles in the chain interior. This behavior is similar for chains in the unbounded 
and the wall-bounded regions. The forces arc smaller for long chains, because the parti- 
cles collectively drag the fluid in the direction of the chain velocity. This mechanism is 
diminished, but not eliminated by the wall presence. 

In an unbounded space, the force distribution in a chain moving in the transverse 
direction y is qualitatively similar to the distribution for the longitudinal motion discussed 
above. In the wall-bounded region the results are, however, considerably different: the 
forces near the center of the chain are much larger than the forces near the chain ends. 
This behavior, clearly seen in figure [TUl for H/2a < 2, stems from the conservation of 
the fluid volume. The chain moving in the transverse direction acts like a piston pushing 
fluid along the space between the walls, thus producing a pressure-driven flow decaying 
on the lengthscale / = 2aN. The pressure increases linearly with the chain length until 
it is large enough to push the fluid back through the gap between the walls and the 
particles. At this point, the pressure becomes independent of N. The pressure produced 
by this mechanism is responsible for the large resistance coefficient of long chains in 
transverse motion between closely spaced walls, as shown in figure El 

For a chain moving in the direction z (normal to the walls) in a system with a small 
value of the wall-particle gap e, the resistance coefficients dominated by the 

lubrication forces between the walls and the individual particles. The coefficients are 
the smallest for the spheres at the chain ends, as seen in figure ^1 unlike for chains in the 
infinite space. This behavior stems from the presence of the geometrical constraints — 
the resistance is smaller where there is more room for the fluid to escape from the gaps 
between the walls and the particles. 

The geometrical-parameter dependence of the torque acting on individual spheres in a 
translating chain is less varied, as illustrated in figure lTTI In all configurations considered, 
we find that the torque on the interior spheres is much smaller than the torque at the 
chain ends. An interesting feature is the sign change of the torque acting on the particle 
i = 2 for the coefficient Cr^. 



8. Conclusions 

Many-body hydrodynamic interactions in suspensions of spherical particles confined 
between two parallel planar walls have been studied here theoretically and numerically . 
Our primary theoretical result is a set of compact expressions for the multipolar matrix 
elements of the Green's integral operator for Stokes flow in the space between the walls. 
The matrix elements are given in the form of lateral Fourier integrals of products of 
several simple matrices. 

Our expressions have been used to develop an algorithm for evaluating many-particle 
hydrodynamic friction and mobility matrices in a wall-bounded suspension. The algo- 
rithm involves solving a set of linear equations for the multipolar moments of the force 
distributions induced on the particles. The resulting friction matrix is corrected for the 
lubrication forces by using a superposition of particle-particle and particle- wall contribu- 
tions. Our algorithm yields highly accurate results — for example, the results presented 
in this paper have been obtained with an accuracy better than 1 %. We note that at 
each truncation of the force-multipole expansion the boundary conditions at the walls 
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are exactly satisfied. This feature is essential for obtaining a proper far-field behavior of 
the friction matrix, including the strong backflow effects seen for rigid arrays of spheres. 

Our numerical algorithm has been used to evaluate the hydrodynamic resistance matrix 
for a single particle, a pair of particles, and a system of many particles confined between 
two planar walls. The problem of hydrodynamic interactions in the two- wall geometry 
involves several characteristic lengthscales: the particle radius a, the wall-to-wall distance 
H, and the lateral distance between the particles p. For p <^ H the interactions between 
particles are similar to those in the infinite space. For p H the crossover occurs to a 
quasi-two-dimensional behavior in the regime p ^ H. 

In the quasi-two-dimensional domain the vertical components of the mutual pair resis- 
tance matrix decay exponentially, and the lateral components behave as 0(/5~^). More- 
over, the sign of the transverse component of the resistance matrix is opposite to the 
sign of this compone nt for a pair of pa rticles in infinite space. As discussed here and 

in our recent paper <)Bhattacharva etal 2005 1 . this behavior can be explained using 

the asymptotic Hele-Shaw (lubrication) form of the far-field flow produced by a moving 
particle. 

The crossover behavior is particularly pronounced for rigid arrays of spheres arranged 
along a line parallel to the walls. In the regime a <^ I <^ H , where I is the chain length, 
the hydrodynamic friction force per particle decreases as {logl)~^ for large I, similar to 
the behavior in the infinite space. In contrast, for I ^ H the longitudinal component of 
the friction tensor (per particle) and the component normal to the walls tend to constant 
values. Moreover, for small particle-wall gaps, the transverse component (normal to the 
chain but parallel to the walls) increases linearly with the chain length before it saturates 
at a value that is much higher than the corresponding value fo r the longitudinal motion. 

As discussed in our recent paper ijBhattacharva et alll2005|) . the standard wall-super- 
position approximation is insufficient for many problems. The resistance matrix in such 
an approximation is composed from two single-wall contributions. In particular, the su- 
perposition approximation yields a wrong sign of the transverse component of the mutual 
pair resistance matrix and a wrong far-field behavior of all components of this matrix. The 
approximation also fails to reproduce the striking increase with the number of particles 
for the transverse resistance coefficient of linear arrays of spheres. 

The numerical efficiency of our method can be substantially improved by combining 
our Cartesian representation of the wall contribution to the Green's matrix with the 
asymptotic far-field expressions for this quantity. The asymptotic expressions which we 
have recently derived can be expressed in terms of multipolar solutions of Laplace's 
equation for two-dimensional pressure field corresponding to the lubrication flow in the 
space between the walls. These expressions do not involve Fourier integrals, and they 
can relatively easily be implemented in numerical algorithms for periodic systems and in 
accelerated PPPM or fast-multipole algorithms. 
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J. B. was supported by NSF grant cts-s0348175 and by Hellman Foundation. 



Appendix A. Spherical basis 

The spherical basis of Stokes flows and the reciprocal basis fields v/f^^ used in 

the present paper are normalized differently than the corresponding basis fields v^^^^^^ 
and W;"^^^^^ introduced bv lCichocki et ali ((l988j) . The relations between these sets of 
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basis functions are as follows: 

_ , . _i -(CFS)/ X + ^ ^ Ar -1 +(CFS). N 
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(Ala) 



wz„<,(r) = A^/<Tn/™rw,^^ ^(r), w+„^(r) = iV^^^ni^rW;^^ ^(r), 



where 



and 



Nio = 1, 
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-{I + Ni2 = l[{l + 1){21 + 1){21 + 



nir. 



An {l + my. 



21 + 1 {I -my. 



1/2 



Below we list the explicit expressions for the angular coefficients V^^^ ( 
ical basis fields in our present normalization, 



(2; + 1)2 

V 



Iml 



and 



A+ - 
lm2 



I 



2(2^+1) 



v+ 



1 



ii-i 



I 



where 



Y,,_i™(f) = arV-'+i-r„i 



(? + l)(2/ + l)(2? + 3) 
V [r'r;„,(f)] 



r/™(r) 



Y/(,„(f) = 7, X Vsy/,„(r) 

are the normalized vector spherical harmonics, as defined bv lEdmondsl l)l96i 



Yi^ii) = n,-;„i(-l)"Pr(cos0)e™'^ 
are the normalized scalar spherical harmonics, and 

ai^[l{2l + l)Y/\ Pi = \{l + l){2l + l)Y/\ 7, = -i[;(/ + l)]i/2. 

Appendix B. Transformation vectors X* and X"^ 

The resistance matrix (I2.1|l is obtained from the solution 

N 



(A1&) 
(A 2) 

(A3) 
for spher- 

(A4a) 

(A4&) 
(A 4c) 
(A 5a) 

(A56) 

(A 5c) 

(A6a) 

(A 66) 

(A 6c) 
. Here 

(A 7) 
(A 8) 



(Bl) 
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of the forcc-multipole equation (|3.8(l by projecting the gencrahzcd friction matrix F = 
M^^ onto the subspaces corresponding to the rigid-body motion of the particle j and the 
total force and torque of the induced-force distribution on particle i. As shown in Ref. I, 
the projection can be expressed in the form 



^^-s = V X(y4 I lma)Fij{lma \ l'm'a')X{l'm'a' \ B), 



Iraa V m' a 



(B2) 



where ^, S = t, r. Here X(A | Ima) and X.{l'm'a' \ B) are the projection vectors defined 
by the equations 



X(t I Ima) = 6ii6cro^^ (m), X(r | Ima) — SnS^i'X.'^im), 



(B3) 



X*(-l) = ( 



—1 




X*(0) = an)' 






V2 



XVl) 



and 



X''(m) = -2iX*(m), to = -1,0,1, 



X{lma \ A) =X.*{A\lma), A = t,r. 



(|7r)V2 



-1 
— i 




(B4) 
(B5) 

(B6) 



Appendix C. Cartesian basis fields 

The Fourier coefficients V^^(z) in the expression H4.1|l for the Cartesian basis fields 
are given by the expressions 



V^o(z) = (327t2)-i/2 i(l - 2kz)k +{1 + 2kz)e, 



-1/2 



-1/2 



and 



Vki(^) = (87t2)-i/2 (k X e,)fc 
V^2(z) = (327r2)-i/2 (ik-e.)r 

V+ (z) = (327r2)-i/2 (ik + e,)fc-i/2, 
V^i(z) = (87t2)-i/2 (k X e,)fc-i/2, 



Vk2(z) = (327t2)-i/2 i(l + 2fcz)k - (1 - 2kz)e, 
where k = h/k. The corresponding pressure fields are 

Pko(l-) = (27t2)-l/2^;jl/2gik-p-fc^^ ^ (27T2)-l/2j^fcl/2gik.p+fcz^ 



and 



Pkl(l") = Pk2(>^) = Pko('') = ?'kl('') = 0- 



(Cla) 
(C1&) 
(Clc) 

(C2a) 
(C2&) 
(C2c) 

(C3) 
(C4) 
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